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Introduction 

Three-component  digital  seismic  refraction  data  were  collected  along  profiles  over  outcrops 
of  the  Troodos  Ophiolite  Complex  (Cyprus)  as  described  by  Dougherty  et  al.  (1992).  Estimates  of 
seismic  P-wave  velocity  for  the  rocks  comprising  the  outcrops  were  obtained  using  three  different 
procedures  applied  to  the  first  arrival  times: 

1 .  Straight-Ray  Procedure :  computation  of  the  reciprocal  of  the  slope  of  a  best-fit  least-squares 
straight  line  through  first  arrival  times  plotted  as  a  function  of  the  distance  along  the  straight 
ray  from  source  to  receiver; 

2.  Hawkins’  Procedure :  computation  of  the  reciprocal  of  the  slope  of  a  best-fit  least-squares 
straight  line  through  corrected  first  arrival  times  plotted  as  a  function  of  source-to-receiver  off¬ 
set  (a  procedure  similar  to  that  outlined  by  Hawkins,  1961,  pp.  809-810); 

3.  Ray-Tracing  Inversion  Procedure :  two-dimensional  ray  tracing  inversion  of  first  arrival  times 
to  determine  a  subsurface  velocity  model  using  the  method  of  Zelt  and  Smith  (1992). 

For  the  second  and  third  procedures,  a  simple  subsurface  velocity  structure  consisting  of  a 

near-surface  low-velocity  zone  overlying  higher  velocity  material  was  assumed  at  each  site.  This 

model  corresponds  to  a  near-surface  weathered  zone  overlying  relatively  unweathered  bedrock. 

The  major  scientific  use  of  the  results  of  the  velocity  analyses  is  the  comparison  of  the  P-wave 

velocities  of  the  unweathered  bedrock  to  the  lithological  differences  between  the  sites.  This  report 

only  addresses  the  velocity  analyses  and  does  not  interpret  the  results  in  lithological  terms. 

Timing  First  Arrivals 

The  first  stage  of  any  interpretation  of  seismic  data  is  the  establishment  of  the  arrival  times  for 
one  or  more  wave  types.  The  first  arrival  time  is  defined  to  be  the  time  (measured  from  the  shot 
initiation  as  time  zero)  of  the  first  obvious  break  of  the  seismic  trace  from  the  background  noise 
level.  This  first  obvious  break  of  the  seismic  trace  corresponds  to  the  arrival  of  the  P-wave  traveling 
between  source  and  receiver  along  the  fastest  available  path.  Thus,  the  first  arrival  time  is  also  com¬ 
monly  called  the  P-wave  travel-time.  Later  arrivals,  which  may  correspond  to  different  wave  types, 
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could  not  be  clearly  discerned  in  the  Troodos  data  because  of  superposition  of  many  different  seis¬ 
mic  signals. 

Timing  of  first  arrivals  was  performed  using  a  MATLAB  routine  after  transferring  the  Bison 
field  records  to  a  DECstation  5000  and  converting  them  to  MATLAB -readable  format  (MATLAB 
is  a  general-purpose  computational  program  for  scientists  and  engineers  and  is  available  from  The 
Math  Works,  Natick,  MA).  The  data  were  not  corrected  for  phase  delays  associated  with  the  seis¬ 
mic  system;  these  delays  were  assumed  to  be  uniform  across  the  spread  length  for  a  given  compo¬ 
nent. 

To  ensure  consistency  in  measurement  of  the  first  arrival  times,  all  the  shot  records  for  a  par¬ 
ticular  site  were  analyzed  by  the  same  person  working  within  an  established  system.  The  compo¬ 
nents  (vertical,  longitudinal,  and  transverse)  for  a  given  shot  were  analyzed  separately. 
Irregularities  in  the  data  in  the  form  of  polarity  reversals  or  the  presence  of  precursors  were  noted. 
After  analysis  of  a  component  for  a  particular  shot  record,  the  first  arrival  times  were  examined  to 
identify  anomalies.  These  anomalies  were  then  checked  for  possible  timing  blunders  and  were  cor¬ 
rected  where  appropriate. 

Straight-Ray  Procedure 

Procedure 

As  a  first  step  in  the  analysis  of  seismic  velocity  at  a  given  Troodos  site,  a  homogeneous  and 
isotropic  elastic  subsurface  was  assumed  so  that  a  single  number  could  be  used  to  describe  the  P- 
wave  velocity.  A  point  source  in  a  material  of  this  type  generates  spherical  wavefronts  and  straight 
rays,  which  in  turn  imply  a  linear  graph  of  first  arrival  time  plotted  as  a  function  of  the  distance 
along  the  straight  ray  from  source  to  receiver.  The  reciprocal  of  the  slope  of  the  linear  graph  gives 
the  seismic  P-wave  velocity  of  the  subsurface. 
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Of  course,  real  seismic  field  data  will  rarely  be  collected  over  a  subsurface  that  closely  approx¬ 
imates  a  homogeneous  and  isotropic  elastic  material  with  a  single  P-wave  velocity.  In  the  case  of 
the  Troodos  seismic  experiments,  it  is  likely  that  there  are  at  least  two  distinct  P-wave  velocities  at 
each  site,  a  relatively  thin  near-surface  low-velocity  layer  representing  weathered  rock,  and  a 
region  of  faster  velocity  representing  competent  shallow  bedrock.  Since  the  near-offset  geophones 
of  the  Troodos  experiments  may  have  recorded  a  P-wave  propagating  through  the  near-surface 
weathered  layer,  the  first  arrival  times  from  near-offset  geophones  were  discarded  to  better  empha¬ 
size  the  P-wave  velocity  of  the  more  competent  shallow  bedrock. 

Computation  of  the  distance  along  the  straight  ray  between  source  and  receiver  was  accom¬ 
plished  using  the  total  station  survey  data  and  the  usual  formula  for  the  distance  between  two  points 
in  space: 

d  =  J(xR  -  x5)2  +  (yR  -  ys)2  +  (zR  -  ZS)2  (EQ1) 

where  d  is  the  distance  between  source  and  receiver,  and  (xs,  ys,  Zs )  and  (xR,  yR,  zR)  are  the  rectan¬ 
gular  coordinates  of  the  source  and  receiver,  respectively. 

After  the  first  arrival-time  data  were  plotted  as  a  function  of  distance  between  source  and 
receiver,  a  best-fit  straight  line  through  the  data  was  computed  using  the  method  of  least  squares. 
The  reciprocal  of  the  slope  of  this  straight  line  constitutes  an  estimate  of  the  P-wave  velocity  of  the 
subsurface.  Scatter  about  the  best-fit  straight  line  indicates  random  timing  errors  and  residual 
shortcomings  in  the  uniform  velocity  assumption.  As  an  example,  Figure  1  shows  the  results  of  a 
velocity  determination  at  Site  1  for  the  longitudinal  first  arrival-time  data  of  Shot  4. 

For  a  given  site,  the  straight-ray  procedure  was  carried  out  separately  for  each  component 
recorded  during  each  shot.  On  center  shots,  velocity  calculations  were  carried  out  in  both  direc¬ 
tions.  Thus,  given  a  full  complement  of  five  shots  at  a  site,  and  assuming  a  single  seismic  profile 
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Figure  1.  Example  of  velocity  determination  using  first  arrival  times  plotted  as  a  function  of  the 
distance  from  source  to  receiver  along  a  straight  ray.  The  graph  is  for  the  longitudinal  component 
of  Shot  4  at  Site  1.  A  velocity  of  1925  m/s  was  calculated  from  the  reciprocal  of  the  slope  of  the 
best-fit  least-squares  straight  line. 


at  the  site,  it  is  possible  to  compute  six  velocities  for  each  geophone  component  and  average  the 
results. 


Results 

A  summary  of  the  P- velocity  analysis  by  the  straight-ray  procedure  is  given  in  the  table  on  the 
next  page.  Sites  are  arranged  according  to  stratigraphic  position  (sites  closest  to  the  paleoseafloor 
given  first)  with  the  geographic  location,  type  of  alteration  zone,  gross  lithology,  and  estimated 
depth  below  the  paleoseafloor  listed  for  each  site  (estimated  depths  below  paleoseafloor  from  Gillis 
and  Sapp,  1997;  depths  for  sites  8-9,  11,  and  12  are  not  known).  At  a  given  site,  the  mean  velocity, 
estimated  standard  deviation  of  the  mean  velocity,  and  velocity  range  are  given  for  each  compo¬ 
nent,  with  the  component  with  the  largest  mean  velocity  denoted  by  boldface  print.  Abbreviations 
used  in  the  table  are  as  follows:  SFWZ-seafloor  weathering  zone;  LTZ-low-temperature  zone;  V- 
vertical  component;  L-longitudinal  component;  T-transverse  component. 
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Site 

Alt 

Zone 

Geographic  Location 

Gross  Lithology 

Depth  Below  Paleoseafloor 

McanPJeUs.d.) 

n 

(m/s) 

1 

SFWZ 

Akaki  Canyon 

a 

2064(116.3) 

4 

2299-1830 

Pillow  Lav2 

a 

2083  (97.2) 

4 

2291-1843 

0  m 

a 

2102  (117.5) 

4 

1817-1317 

3 

SFWZ 

Akaki  Canyon 

a 

2395  (84.0) 

8 

2821-2085 

Pillow  Lava 

L 

2326  (63.9) 

8 

2696-2093 

0  m 

T 

2205  (66.3) 

8 

2425-1906 

2 

SFWZ 

Akaki  Canyon 

a 

1537  (146.3) 

3 

1814-1317 

Pillow  Breccia 

a 

1534  (128.9) 

3 

1775-1274 

50  m 

a 

1468  (145.0) 

3 

1784-1170 

6 

LTZ 

Margi 

a 

2813  (324.1) 

6 

4391-2263 

Differentiated  Massive  Flow 

a 

2488  (62.8) 

6 

2781-2385 

50-100  m 

a 

2369(121.2) 

6 

2930-2157 

MM 

LTZ 

Akaki  Canyon 

a 

2471  (276.0) 

5 

3995-1878 

Pillow  Lava 

a 

2301  (366.2) 

5 

3710-1571 

1  1 

200  m 

a 

2362  (387.8) 

5 

3829-1527 

10 

LTZ 

Mathiati  Mine 

a 

2208  (44.9) 

4 

2309-2111 

Pillow  Lava 

2162  (44.2) 

4 

2294-2106 

300-400  m 

T 

2008  (86.3) 

4 

2136-1755 

H 

LTZ 

Analiondas 

a 

3469  (267.9) 

8 

4553-2666 

SI 

Pillow  Lava 

a 

8 

1  1 

500  m 

T 

3154  (205.2) 

8 

4225-2607 

5 

LTZ 

Eucalyptus  Canyon 

a 

2813  (626.8) 

3 

4065-2120 

Massive  Flow 

L 

2580  (369.2) 

3 

3318-2191 

500  m 

a 

2642  (433.5) 

3 

3503-2127 

8-9 

LTZ 

Kampia  Mine 

a 

2726  (229.8) 

8 

4298-2265 

Intercalated  Hyaloclastites,  Pillow 

L 

2551  (86.8) 

8 

2905-2210 

Lavas,  Sheet  Flows,  Sheeted  Dikes 
??  m 

T 

! 

2422  (91.7) 

8 

2908-1956 

11 

LTZ 

Delikipo  Area 

a 

3019  (92.0) 

6 

3191-2584 

Pillow  Lavas  Intermixed  with  Thin 

a 

3118  (64.7) 

6 

3380-2905 

Sheet  Flows 
??  m 

a 

2689  (124.9) 

6 

3047-2163 

12 

LTZ 

Peristerona  River 

m 

Straight-ray  computations  not  carried  out  at  Site  12. 

Series  of  Sheet  Flows  Underlain  by 

■ 

Single-component  geophones  were  deployed  in  soil 

Pillow  Lavas 

T 

overlying  bedrock. 

??  m 

Examination  of  the  standard  deviations  in  the  preceding  table  indicates  that  at  a  particular  site, 
the  differences  between  the  mean  velocity  measured  on  the  vertical  component  and  the  mean  veloc¬ 
ities  measured  on  the  horizontal  components  are  not  statistically  significant.  However,  if  one 
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examines  the  collection  of  all  sites,  the  average  vertical  velocity  is  larger  than  the  average  horizon¬ 
tal  velocities  (longitudinal  and  transverse)  at  8  of  10  sites  (the  only  exceptions  are  Sites  1  and  11). 
This  observation  suggests  statistical  testing  of  the  null  hypothesis  that  there  exists  no  bias  towards 
measuring  the  maximum  velocity  on  the  vertical  component.  The  formal  statement  of  the  null 
hypothesis  is: 

H0:  p  =  1/3  where  p  is  the  probability  of  measuring  the  maximum  velocity 
on  the  vertical  component, 

In  other  words,  H0  states  that  in  any  given  velocity  determination,  the  maximum  velocity  is  just  as 
likely  to  be  measured  on  any  of  the  three  components,  so  that  the  probability  of  measuring  the  max¬ 
imum  on  the  vertical  component  is  one  third.  The  appropriate  test  statistic  is  given  by  Pollard 
(1977,  pp.  137-139): 


(r-np0) 

Jnp0{\-p0) 


(EQ  2) 


where  T  is  approximately  distributed  normally  with  zero  mean  and  unit  variance,  p0  =  1/3  (the 
probability  of  measuring  the  maximum  velocity  on  the  vertical  component  if  there  is  no  bias),  r  is 
the  number  of  trials  where  the  maximum  velocity  is  measured  on  the  vertical  component  ( r  =  28  in 
this  case),  and  n  is  the  total  number  of  trials  for  the  entire  Troodos  data  set  (n  =  55  which  is  the 
number  of  shots  recorded  for  which  straight-ray  velocities  were  determined  using  all  three  compo¬ 
nents).  The  value  of  T is  2.765  and  standard  tables  for  the  unit  normal  distribution  indicate  that  the 
probability  of  a  larger  value  is  less  than  0.003.  Thus,  unless  a  rare  outcome  has  occurred,  the  null 
hypothesis  is  rejected  at  a  very  high  level  of  significance. 


The  outcome  of  the  hypothesis  test  suggests  that  there  is  an  overall  tendency  to  time  the  first 
arrival  on  the  vertical  component  earlier  than  on  the  horizontal  components  of  the  same  geophone. 
Although  this  result  was  not  confirmed  by  a  rigorous  check  of  the  entire  set  of  first  arrival  times,  it 
is  consistent  with  the  genera]  pattern  noted  during  the  timing  stage  of  the  data  analysis.  The  sim- 
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plest  explanation  is  a  difference  in  system  response  between  the  horizontal  and  vertical  compo¬ 
nents  (i.e.,  a  difference  in  phase  delay  between  the  vertical  and  horizontal  components  with  the 
vertical  component  early).  However,  it  is  also  possible  that  the  particle  velocity  at  the  geophone 
during  the  first  arrival  may  be  oriented  more  vertically  than  horizontally,  thereby  giving  a  more 
impulsive  signature  on  the  vertical  component  that  tends  to  be  timed  earlier. 

The  particle  motion  of  first  arrivals  for  the  Troodos  seismic  experiments  supports  the  conclu¬ 
sion  that  the  first-arriving  seismic  energy  involves  particle  velocities  that  are  more  vertical  than 
horizontal.  Specifically,  633  of  912  particle  motion  diagrams  (69.4%)  indicate  an  initial  angle  of 
particle  velocity  of  45°  or  greater  with  the  horizontal.  Is  the  predominance  (nearly  70%)  of  more 
vertically  oriented  particle  velocity  consistent  with  the  basic  assumption  of  P-waves  propagating 
along  straight  rays?  This  question  is  not  addressed  in  this  report  because  it  is  complicated  by  the 
topography  of  the  Troodos  sites  (see  Appendix)  and  the  partitioning  of  seismic  energy  at  the  free 
surface.  Furthermore,  rather  than  assuming  straight  rays,  it  may  be  more  correct  to  view  first-arriv¬ 
ing  seismic  energy  as  following  the  path  of  a  P-head  wave  generated  at  the  base  of  the  weathered 
layer  and  then  returning  to  the  surface  at  a  steep  angle.  It  is  this  latter  possibility  that  is  explored 
in  the  next  sections  (Hawkins’  procedure  and  ray-tracing  inversion  procedure). 

Hawkins’  Procedure 

Procedure 

The  term  Hawkins'  procedure  refers  here  to  a  procedure  for  velocity  determination  that  is 
associated  with  the  reciprocal  method  of  seismic  refraction  interpretation  (Hawkins,  1961).  The 


Figure  2.  Subsurface  model  and  rays  used  by  Hawkins  (1961)  to  develop  the  reciprocal  method. 
Sources  A  and  B  are  at  the  surface  with  receiver  at  G .  Point  P  is  the  midpoint  of  the  base  of  an 
isosceles  triangle  formed  by  the  upgoing  rays  to  G. 
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subsurface  is  assumed  to  consist  of  a  near-surface  low-velocity  material  (Va)  that  overlies  material 
of  higher  velocity  (Vj)  so  that  a  head  wave  can  be  generated  at  the  interface  (Figure  2).  Consider 
the  reversed  refraction  profile  from  source  A  to  source  B,  and  visualize  three  head-wave  rays: 


•  the  ray  from  A  to  receiver  G, 

•  the  ray  from  B  to  receiver  G,  and 

•  the  ray  from  source  A  to  source  B. 

Now  let  tAG  be  the  travel  time  from  A  to  G,  let  tBG  be  the  travel  time  from  B  to  G,  and  let  tAB 
be  the  travel  time  from  A  to  B.  Using  this  notation,  a  new  quantity,  the  corrected  travel  time 
( t'AG ),  may  be  defined: 


t’ 


tAG~tBG  +  fAB 


(EQ  3) 


MG  -  2  ' 

The  usefulness  of  t' AG  is  that  it  may  be  shown  to  equal  the  travel  time  from  A  to  P  along  the 

partial  head-wave  ray  path  ALP  (Hawkins}  1961): 


.  _  AL  LP 

AG  y0  vl 


(EQ  4) 


To  a  first  approximation,  the  distance  along  LP  may  be  represented  by  the  offset  x  measured 
along  the  surface  between  A  and  G: 


,  _  AL  x 

AG~yo  +  y^ 


(EQ  5) 


Equation  (5)  indicates  that  t'AG  is  a  linear  function  of  the  offset  x  with  slope  1 IV j  and  intercept 
AL/V0.  Thus,  given  the  travel  times  required  by  equation  (3)  for  a  number  of  receiver  positions 
between  sources  A  and  B,  a  plot  of  ^  AG  as  a  function  of  offset  x  will  yield  a  straight  line  whose 
slope  can  be  inverted  to  give  the  velocity  Vj  (i.e.,  the  velocity  of  the  material  beneath  the  near¬ 
surface  low-velocity  layer). 


The  application  of  Hawkins’  procedure  to  the  Troodos  seismic  data  requires  that  Vc>  be  identi¬ 
fied  with  the  relatively  thin  near-surface  low- velocity  layer  representing  weathered  rock,  and  Vj  be 


identified  with  the  underlying  region  of  faster  velocity  representing  competent  shallow  bedrock. 
Since  each  site  had  multiple  pairs  of  shots  to  represent  sources  A  and  B,  it  was  possible  to  repeat 
the  t'AG  analysis  several  times  for  a  given  site  and  then  compute  a  mean  P-wave  velocity  for  that 
site.  Because  the  first  arrival  times  tend  to  be  earlier  on  the  vertical  component,  the  Hawkins’  pro¬ 
cedure  was  only  applied  to  data  recorded  on  the  vertical  components. 

Results 

A  summary  of  the  P-velocity  analysis  by  the  Hawkins’  procedure  is  given  in  the  table  below. 
As  in  the  previous  velocity  table,  sites  are  arranged  according  to  stratigraphic  position  (sites  closest 
to  the  paleoseafloor  given  first)  with  the  geographic  location,  type  of  alteration  zone,  gross  lithol¬ 
ogy,  and  estimated  depth  below  the  paleoseafloor  listed  for  each  site  (estimated  depths  below  sea¬ 
floor  from  Gillis  and  Sapp,  1997;  depths  for  sites  8-9,  11,  and  12  are  not  known).  At  a  given  site, 
the  mean  velocity  and  estimated  standard  deviation  of  the  mean  velocity  are  given  for  both  the 
Hawkins’  procedure  and  the  straight-ray  procedure;  results  from  vertical  component  data  only  are 
shown.  The  larger  of  the  two  mean  velocities  at  a  given  site  is  given  in  boldface  print.  Abbrevia¬ 
tions:  SFWZ-seafloor  weathering  zone;  LTZ-low-temperature  zone. 


Site 

Alt 

■  Zone 

Site  Location  \  ~ 

,  Gross  Lithology  ■  '  Vv 
j  p-  Depth  Below  Paleoseafloor 

Hawkins’  Procedure 
Mean  P-Vel  ($.  d  )  Sfl 
No.  Meas. 

r  ••••.-. (m/s)  *  - 

Straight-Ray  Proc.  (Vert) 
Mean  Vel  (s.  d.) 
No.Meas. 

(ni/s)  y::’v 

1 

SFWZ 

Akaki  Canyon 

Pillow  Lava 

0  m 

2090  (92.5) 

4 

2064(116.3) 

4 

3 

SFWZ 

Akaki  Canyon 

Pillow  Lava 

Om 

2420  (50.7) 

8 

2395  (84;0) 

8 

2 

SFWZ 

Akaki  Canyon 

Pillow  Breccia 

50  m 

1640  (140.2) 

2 

1537(146.3) 

3 

6 

LTZ 

Margi 

Differentiated  Massive  Flow 
50-100  m 

2610(63.2) 

6 

2813  (324.1) 

6 

9 


!PH 

ii 

(S'?i 

■■m om» 

Site  Location 

■  '3fl 

: 

Ml  '/M 

iSf-M 

mm 

1  Strai^ht-R! *y  Proe  (Veit) 

r'>„. 

If 

DepthBelow  Paleoseafloor 

'  j  „ 

(m/ )  it 

■ 

LTZ 

Akaki  Canyon 

Pillow  Lava 

200  m 

2780  (115.5) 

2 

2471  (276.0) 

5 

10 

LTZ 

Mathiati  Mine 

Pillow  Lava 

300-400  m 

2180  (67.6) 

4 

2208  (44.9) 

4 

■ 

LTZ 

Analiondas 

Pillow  Lava 

500  m 

3120  (56.5) 

8 

3469  (267.9) 

8 

5 

LTZ 

Eucalyptus  Canyon 

Massive  Flow 

500  m 

2480(134.9) 

2 

2813  (626.8) 

3 

8-9 

LTZ 

Kampia  Mine 

Intercalated  Hyaloclastites,  Pillow 
Lavas,  Sheet  Flows,  and  Sheeted 
Dikes 
??  m 

2510(49.2) 

8 

2726  (229.8) 

8 

11 

LTZ 

i 

Delikipo  Area 

Pillow  Lavas  Intermixed  with  Thin 
Sheet  Flows 
??  m 

3030  (83.8) 

8 

3019  (92.0) 

6 

12 

LTZ 

Peristerona  River 

Series  of  Sheet  Flows  Underlain 
by  Pillow  Lavas 
??  m 

i 

2240  (62.9) 

4 

2790  (83.5) 

4 

Straight-ray  computations 
not  carried  out  at  Site  12. 

!  Single-component  geophones 
were  deployed  in  soil  overly¬ 
ing  bedrock. 

Comparison  of  results  from  the  straight-ray  and  Hawkins’  procedures  indicates  reasonable 
consistency  in  P-wave  velocity  determinations.  This  statement  is  strengthened  if  one  considers  the 
estimated  standard  deviations  of  the  mean  velocities.  There  is  no  overall  tendency  for  mean  veloc¬ 
ity  estimates  to  be  higher  for  the  Hawkins’  procedure  (i.e.,  the  boldface  velocities  occur  5  times  in 
the  Hawkins’  column  and  5  times  in  the  straight-ray  column). 


Ray-Tracing  Inversion  Procedure 

Ray  tracing  may  be  combined  with  modem  inversion  theory  to  provide  a  powerful  velocity 
estimation  technique.  In  this  technique,  observed  first  arrival  times  are  compared  with  theoretical 
times  computed  according  to  geometrical  ray  theory  applied  to  a  realistic  subsurface  velocity 
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model.  The  initial  model  is  systematically  modified  until  the  fit  between  the  observed  and  theoret¬ 
ical  travel  times  meets  some  predetermined  criteria. 


Ray-tracing  inversion  of  the  Troodos  seismic  data  was  accomplished  using  RAYINVR,  a  pub¬ 
lic  domain  program  developed  at  the  University  of  Utah  (Zelt  and  Smith,  1992)  and  based  on  a  fast 
raytracing  algorith  developed  at  the  University  of  British  Columbia  (Zelt  and  Ellis,  1988).  A  fun¬ 
damental  feature  of  RAYINVR  is  its  two-dimensional,  layered,  variable-sized  block  parameteriza¬ 
tion  which  allows  for  both  efficient  ray  tracing  and  a  minimum  number  of  independent  model 
parameters.  Although  RAYINVR  was  developed  for  use  on  crustal  seismic  data  (scale  of  profiles 
measured  in  hundreds  of  km),  the  program  is  very  flexible  and  can  be  used  effectively  for  outcrop 
seismic  data  (scale  of  profiles  less  than  100  m).  The  theory  and  methodology  on  which  RAYINVR 
is  based  are  described  by  Zelt  and  Smith  (1992). 

Methodology  for  Ray  Tracing 

The  model  is  parameterized  as  a  sequence  of  layers,  and  each  layer  is  comprised  of  a  series  of 
adjacent  trapezoids  with  vertical  sides.  Velocities  may  be  continuous  or  discontinuous  across  layer 
boundaries  but  are  continuous  laterally  within  a  given  layer.  Layers  may  be  reduced  to  zero  thick¬ 
ness  so  that  pinchouts  and  bodies  of  finite  lateral  dimension  can  be  modeled.  In  order  to  minimize 
scattering  and  focusing  of  rays  associated  with  refraction  at  layer  boundaries  (which  consist  of 
straight  line  segments),  RAYINVR  includes  an  option  that  simulates  smooth  layer  boundaries. 


Ray  tracing  within  a  layer  is  carried  out  by  solving  one  of  the  following  sets  of  equations  using 
a  Runge-Kutta  method: 


=  cot0 
ax 

^  =  COt0 

dz 


dd  _  (vz-V  cote) 

dx  v 

dQ  _  (vz  ■  tan0  -  vx) 
dz  v 


(for  near  horizontal  ray  paths) 
(for  near  vertical  ray  paths) 


(EQ  6) 

(EQ7) 
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where  x  and  z  are  the  horizontal  and  vertical  rectangular  coordinates,  respectively,  0  is  the  angle 
between  the  tangent  to  the  ray  and  the  z-axis,  v  is  the  wave  velocity  at  point  (x,  z),  and  the  partial 
derivatives  of  v  with  respect  to  x  and  z  are  given  by  vx  and  vv  respectively.  The  initial  conditions 
require  that  the  starting  point  (x0,  z0)  and  take-off  angle  (0O)  of  the  ray  be  specified.  Ray  step 
lengths  (increments  in  either  the  x  or  z  directions)  are  automatically  determined  to  improve  effi¬ 
ciency  without  sacrificing  accuracy.  Snell’s  law  is  applied  at  the  intersection  of  a  ray  with  a  layer 
boundary.  Rays  may  be  traced  that  correspond  to  turning  waves,  reflected  waves,  and  head  waves. 
Methodology  for  Inversion 

The  formulation  of  P-wave  travel-time  inversion  as  implemented  by  RAYINVR  begins  with 
the  following  assumptions: 

•  the  propagation  of  P- waves  through  the  subsurface  is  accurately  represented  by  zero-order 
asymptotic  ray  theory  in  a  two-dimensional  velocity  model; 

•  the  two-dimensional  velocity  model  can  be  adequately  described  by  M  model  parameters, 
where  M  is  a  finite  number. 

Suppose  we  have  Admeasured  P-wave  travel-times  from  a  given  seismic  experiment.  If  tt  represents 
the  theoretical  travel  time  for  source-receiver  pair  i,  then  the  assumptions  listed  above  imply  that  ti 
is  a  function  of  Mmodel  parameters  (denoted  in  this  report  by  mj,  m2, ...,  mM),  and  that  tt  will  equal 
the  observed  travel  time  when  the  model  parameters  take  on  their  “true”  values  (the  parameters  that 
represent  the  actual  subsurface).  Now  let  Hi,  r|2>  •••>  H m  be  the  “true”  model  parameters,  and  let  fj.j, 
|12,  •••,  he  the  “trial”  model  parameters  (any  parameters  other  than  the  true  parameters).  Using 
these  definitions,  the  observed  travel  time  at  station  i  can  be  expressed  as  a  Taylor  series  expansion 
of  t{  about  the  trial  model  parameters: 

M  ^ 

^1,1)2.  M-2 . M'A#)  +  X  +0  (EQ9) 

J  =  1  7  ”h  =  no¬ 

where  O  represents  higher  order  terms  and  i  =  1,2, N.  If  the  trial  model  is  sufficiently  close  to 
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the  true  model,  then  the  higher  order  terms  may  be  neglected  and  equation  (9)  can  be  rewritten  in 
linearized  form  as  follows: 


M 

At i  =  ^  Am 
7=1 


■idmj 


mk  =  a* 


i  =  1,2,  ...,N 


(EQ  10) 


where  Any  =  (tj  •  -  (iiy)  is  the  parameter  correction  vector  and  At,-  is  the  travel-time  residual  given  by: 


At,-  =  tj-Crij,  rj2, t]M)  -  \i2, \iM)  i  =  1,2,  ...,N 


(EQ  11) 


The  matrix  version  of  equation  (10)  is: 

At  =  A  Am  (EQ  12) 

where  A  is  the  N  x  M  matrix  of  partial  derivatives.  The  model  parameters  used  by  RAYINVR  are 
the  z-coordinates  of  special  points  on  each  surface  between  layers  (i.e.,  boundary  nodes)  and  the 
P-wave  velocities  at  special  points  at  the  top  and  bottom  of  each  layer  (i.e.,  velocity  nodes). 


The  solution  of  equation  (12)  for  Am  is  a  discrete  inverse  problem  (see  Menke,  1989,  for  a 
popular  monograph  on  discrete  inverse  theory).  A  basic  feature  of  this  problem  is  that  it  is  common 
for  some  model  parameters  to  be  overdetermined  and  some  to  be  underdetermined.  Under  these 
conditions,  the  methods  of  discrete  inverse  theory  define  the  solution  vector  (denoted  by  Ams)  to 
be  the  parameter  correction  vector  that  simultaneously  minimizes  the  norm  of  the  prediction  error 

||At-AAms||  (EQ  13) 


plus  the  norm  of  the  parameter  correction  vector  (i.e.,  the  solution  length): 

T 

Ant  Arn  . 
s  s 


(EQ  14) 


If  the  L2  norm  is  selected,  then  simultaneous  minimization  of  (13)  plus  (14)  yields  the  damped 
least-squares  solution: 


where  D  is  the  damping  factor  (determines  the  trade-off  between  minimizing  prediction  error  and 
solution  length),  Cf  and  Cm  are  the  estimated  data  and  model  covariance  matrices,  respectively, 

C  t  =  diag{oj}  C  m  =  diag{o2j},  (EQ 16) 

and  where  a,  is  the  estimated  uncertainty  of  /,•  (considering  equipment  limitations  and  human  error 
in  timing  first  arrivals)  and  Gj  is  an  a  priori  estimate  of  uncertainty  of  parameter  my  The  inverses 
of  Ct  and  Cm  in  (15)  serve  as  weighting  matrices  so  that  some  observations  and  model  parameters 
can  be  emphasized.  The  text  by  Menke  (1989)  and  the  paper  by  Zelt  and  Smith  (1992)  may  be 
consulted  for  details  on  the  damped  least  squares  method,  including  the  question  of  the  trade-off 
between  resolution  and  variance. 

Ray-tracing  inversion  of  travel  times  is  nonlinear  in  the  sense  that  the  higher  order  terms  in 
equation  (9)  cannot  be  neglected  in  general.  The  standard  practice  is  to  choose  a  reasonable  set  of 
initial  trial  parameters  and  then  solve  the  nonlinear  inversion  by  successive  linear  iterations.  A 
formulation  of  the  basic  inversion  procedure  is  as  follows: 

1.  Estimate  the  uncertainties  <7,  and  o;  and  create  the  covariance  matrices  (equation  16). 

2.  Choose  initial  trial  model  parameters  p1?  p2,  •••>  V^m- 

3.  Compute  At  using  equation  (1 1). 

4.  Compute  the  matrix  of  partial  derivatives  A  evaluated  at  pls  p2, ...,  p^  (the  method  used  in 
RAYINVR  is  described  by  Zelt  and  Smith,  1992). 

5.  Compute  Ams  according  to  damped  least-squares  (equation  15)  with  a  damping  parameter  typ¬ 
ically  chosen  between  1-10. 

6.  Add  Ams  to  the  initial  trial  model  parameters  to  obtain  a  new  trial  model. 

7.  Repeat  steps  3-6  until  the  root-mean-square  (RMS)  travel-time  residual  is  within  specified  lim¬ 
its. 
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Initial  RAYINVR  Trial  Model 

The  upper  or  ground  surface  of  the  initial  RAYINVR  trial  model  for  a  given  Troodos  site  was 
defined  by  total  station  topographic  survey  data  (Dougherty  et  al.,  1992).  One  additional  surface 
was  defined  at  a  constant  depth  beneath  the  ground  surface.  This  additional  surface  was  intended 
to  represent  a  weathered  layer  overlying  a  layer  of  competent  bedrock.  Results  from  the  straight- 
ray  and  Hawkins’  procedures  (see  preceding  tables)  were  used  to  determine  initial  trial  velocity 
parameters.  Minor  adjustments  were  made  until  a  reasonable  fit  of  observed  to  calculated  travel 
times  was  achieved.  At  this  point,  the  initial  trial  model  was  declared  ready  for  an  inversion  by  the 
damped  least  squares  method  outlined  above.  Figure  3  shows  an  example  of  a  RAYINVR  initial 
trial  model  for  Site  8  at  Kampia  mine. 

DISTANCE  (m) 


DISTANCE  (m) 


Figure  3.  Initial  trial  model  for  Site  8,  located  on  the  north  face  of  Kampia  mine.  The  upper 
portion  of  the  diagram  illustrates  the  initial  trial  model,  with  ray  paths  for  direct  waves  and 
head  waves  also  shown.  Dashed  lines  indicate  the  vertical  boundaries  of  the  model  blocks. 
Initial  trial  velocity  and  depth  values  are  entered  for  each  corner  of  these  blocks.  The  lower 
plot  shows  the  observed  first  arrival  times  as  crosses  with  the  calculated  first  arrival  times  as 
solid  lines.  The  first  layer  is  uniformly  2-m  thick  for  the  entire  line,  with  a  uniform  velocity  of 
1.25  km/s.  Velocity  of  the  second  layer,  2.50  km/s,  was  taken  directly  from  results  of  the 
Hawkins’  procedure.  This  initial  trial  model  fits  the  data  reasonably  well,  with  calculated 
times  being  slightly  greater  than  real  times.  Only  the  ray  paths  for  one  shot  are  shown  here, 
whereas  in  the  inversion  all  shots  at  Site  8  are  considered. 
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Inversion 

The  RAYINVR  inversion  of  P-wave  travel  time  data  at  a  particular  site  involve  changes  in 
depth  of  the  boundary  nodes,  changes  in  P-velocity  at  the  velocity  nodes,  or  simultaneous  changes 
in  both  (Figures  4-6).  Changes  must  not  violate  known  geology  or  the  measured  topography; 


DISTANCE  (m) 

0  18  36  54  72  90 


Figure  4.  Depth  variant  inversion  of  the  initial  model  shown  in  Figure  3  for  Site  8.  All  velocities 
are  kept  constant.  Note  that  there  are  now  two  distinct  depressions  in  the  upper  layer. 

unwanted  changes  can  be  controlled  by  a  RAYINVR  option  that  permits  fixed  values  at  boundary 
and  velocity  nodes.  Maximum  changes  in  depth  and  velocity  per  iteration  were  set  at  0.05  m  and 
0.05  km/s,  respectively.  As  RAYINVR  converges  on  a  solution,  some  of  the  observed  first  arrival 
times  will  be  discarded  by  the  program  if  they  cannot  be  reasonably  incorporated  into  the  inversion. 

The  final  inversions  indicate  that  a  low  velocity  layer  is  present  at  the  surface  at  all  the  Troodos 
sites.  Models  without  a  surficial  low-velocity  layer  failed  to  match  the  arrival  times  satisfactorily 
and  often  introduced  unreasonable  velocities  or  changes  in  layer  thickness  that  could  not  be  recon¬ 
ciled  with  the  outcrop  geology.  The  low  velocity  layer  probably  represents  the  surficial  weathering 
zone  noted  by  Gillis  and  Sapp  (1997)  at  the  Troodos  sites.  At  two  sites  (6  and  12),  attempts  to  use 
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a  two-layer  model  repeatedly  failed,  but  once  a  third  layer  was  introduced  the  inversion  produced 
a  very  close  match  to  the  observed  travel  times.  In  each  case  the  introduction  of  a  third  layer  makes 
geological  sense  (at  Site  6,  the  middle  layer  appears  to  represent  the  colonnade  at  the  site,  overlain 
by  less  competent  rock;  at  Site  12,  the  layers  represent  overburden,  weathered  rock,  and  more  com¬ 
petent  rock).  All  other  sites  were  modeled  using  two  velocity  layers. 

There  is  an  inherent  weakness  in  the  final  RAYENVR  models.  At  each  site,  the  velocity  of  the 
first  layer  is  based  on  only  a  few  direct  arrivals.  Trade-offs  exist  between  the  thickness  of  the  first 
layer  and  its  velocity.  Slow  velocities  will  result  in  a  thinner  first  layer,  and  faster  velocities  in  a 
thicker  first  layer.  The  mapped  outcrop  geology  is  the  best  check  on  whether  the  first  layer  prop¬ 
erties  are  reasonable.  Even  so,  in  many  areas  the  geologic  information  has  limitations,  and  the 
thickness  and  velocity  of  the  first  layer  is  open  to  question  at  many  sites. 
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DISTANCE  (m) 
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Figure  6.  Combination  inversion  with  variations  in  depth  and  velocity  of  the  initial  model  shown  in 
Figure  3  for  Site  8.  For  clarity,  only  rays  from  one  shot  are  shown.  The  changes  in  the  thickness  of 
the  low  velocity  layer  are  as  in  the  depth  inversion  (Figure  4),  while  the  velocity  changes  in  the  sec¬ 
ond  layer  are  more  subdued  than  in  the  velocity  inversion.  The  depth  changes  are  reasonable  when 
the  geology  of  the  site  is  taken  into  consideration,  so  this  inversion  can  be  viewed  as  an  acceptable 
interpretation  of  the  seismic  data  for  this  site. 


Final  Models 

The  final  RAYINVR  models  for  the  Troodos  sites  are  shown  in  the  Appendix.  All  length 
scales  are  in  meters  (m),  and  all  time  scales  are  in  milliseconds  (ms).  Velocities  at  different  points 
in  the  final  models  are  given  in  kilometers  per  second  (km/s). 
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Appendix 

Final  RAYINVR  Models 


Site  2  RAYINVR 


Site  3c  RAYINVR 


Site  5  RAYINVR 


Site  61  RAYINVR 


Site  6u  RAYINVR 


Site  8  RAYINVR 


Site  9  RAYINVR 


Site  12u  RAYINVR 


